﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using DotSpatial.Controls;
using DotSpatial.Data;
using DotSpatial.Topology;

using DotSpatial.Projections;
namespace GeostatisticalTool.Classes
{
    public class Clip
    {
        private IFeatureSet poly;
        private VectorGrid vector;
        private double cellsize;
        private double[,] data;
        private Extent currentExtent;
        private bool test = false;
        private FeatureSet li;
        public Clip(IFeatureSet poly)
        {
            this.poly = poly;

        }

        public Coordinate GetXYCell(IRaster raster, int col, int row)
        {
            double xmin = raster.Xllcenter;
            double ymin = raster.Yllcenter;
            double cell = raster.CellWidth;
            int numColumns = raster.NumColumns;
            int numRows = raster.NumRows;

            double y = ((numRows - row) * cell) + ymin;
            double x = (col * cell) + xmin;
            return new Coordinate(x, y);
        }

        public IRaster Surface(IRaster raster)
        {



            int nr = raster.NumRows;
            int nc = raster.NumColumns;

            for (int row = 0; row < nr; row++)
            {
                for (int col = 0; col < nc; col++)
                {

                    raster.Value[row, col] = data[col, row];
                }
            }
            raster.Save();
            raster.Close();
            return raster;
        }



        /// <summary>
        /// This function creates the HDR of Gridfile
        /// </summary>
        /// <param name="inExtents"> Extension of grid</param>
        /// <param name="cellSize">Size cell of the grid</param>
        /// <param name="projection">Projection (the same that shapefile)</param>
        /// <param name="noDataValue">No value definition</param>
        /// <param name="outGridPath">Path of the output</param>
        /// <param name="outGrid">Name of the output grid</param>
        public static void CreateGridFromExtents(
            Extent inExtents, double cellSize, ProjectionInfo projection, double noDataValue, string outGridPath, out IRaster outGrid)
        {

            double height = Math.Abs(inExtents.MaxY - inExtents.MinY);
            double width = Math.Abs(inExtents.MaxX - inExtents.MinX);
            int numberRows = Convert.ToInt32(Math.Ceiling(height / cellSize)) + 1;
            int numberCols = Convert.ToInt32(Math.Ceiling(width / cellSize)) + 1;

            // outGrid = Raster.CreateRaster(@outGridPath, null, demRaster.NumColumns, demRaster.NumRows, 1, demRaster.DataType, rasterOptions);
            outGrid = Raster.CreateRaster(@outGridPath, null, numberCols, numberRows, 1, typeof(float), new String[] { });
            outGrid.NoDataValue = noDataValue;
            outGrid.Projection = projection;
            outGrid.CellHeight = cellSize;
            outGrid.CellWidth = cellSize;
            //if (inExtents.MinX < 0)
            //    outGrid.Xllcenter = inExtents.MinX + (cellSize / 2.0);
            //else
            outGrid.Xllcenter = inExtents.MinX;// -(cellSize / 2.0);

            //if (inExtents.MinY < 0)
            //    outGrid.Yllcenter = inExtents.MinY + (cellSize / 2.0);
            //else
            outGrid.Yllcenter = inExtents.MinY;// -(cellSize / 2.0);
        }


        public IRaster GetClip(Extent e, double cellsize, string output)
        {
            currentExtent = e;
            if (test)
            {
                li = new FeatureSet(FeatureType.Line);
                li.Projection = poly.Projection;
                li.DataTable.Columns.Add(new System.Data.DataColumn("id"));
            }

            this.cellsize = cellsize;


            AreaInterest a = new AreaInterest();

            a.Azimut = 0;
            a.CellSizeX = cellsize;
            a.CellSizeY = cellsize;
            a.MinX = e.MinX;
            a.MinY = e.MinY;
            a.MaxX = e.MaxX;
            a.MaxY = e.MaxY;

            IRaster raster;
            // this.np=this.Search.MaxNumPoints;
            CreateGridFromExtents(e, cellsize, poly.Projection, -9999, output, out raster);
            a.NumColumns = raster.NumColumns;
            a.NumRows = raster.NumRows;
            vector = new VectorGrid(a, null);
            //  vector.UpdateAreaInterestByCellSize();
            data = new double[vector.Get_area().NumColumns, vector.Get_area().NumRows];

            foreach (Feature fea in poly.Features)
            {
                InterceptLine(fea);
            }
            if (test)
                li.SaveAs("c:\\temp\\Lakelines.shp", true);

            return Surface(raster);


        }
        private void InterceptLine(Feature fea)
        {


            Extent exf = fea.Envelope.ToExtent();

            if (currentExtent.Contains(exf))
            {
                #region   exteny
                int endRow = vector.GetCell(exf.MinX, exf.MinY)[1];
                int iniRow = vector.GetCell(exf.MaxX, exf.MaxY)[1];

                if (endRow == -1 || iniRow == -1) return;


                for (int it = iniRow; it <= endRow; it++)
                {
                    Coordinate ini = new Coordinate();
                    Coordinate end = new Coordinate();
                    double[] i = vector.GetCoordenatesByCell(1, it);
                    double[] e = vector.GetCoordenatesByCell(vector.Get_area().NumColumns, it);
                    ini.X = i[0] - cellsize;
                    ini.Y = i[1];
                    end.X = e[0] + cellsize;
                    end.Y = e[1];
                    LinearRing line = new LinearRing(new List<Coordinate>() { ini, end });

                    if (test)
                    {
                        IFeature f = li.AddFeature((IBasicGeometry)line);
                        f.DataRow["id"] = it.ToString();
                    }

                    if (fea.Intersects(line))
                    {
                        IFeature rfea = fea.Intersection(line);
                        IList<Coordinate> cor = rfea.Coordinates;
                        if (cor.Count == 1) break;
                        Filldata(rfea, cor);


                    }
                    //else if (fea.Contains(line))
                    //{

                    //    IList<Coordinate> cor = fea.Coordinates;
                    //    Filldata(fea, cor);

                    //}


                }


                #endregion
            }
            else
            {
                for (int it = 0; it < vector.Get_area().NumRows; it++)
                {
                    Coordinate ini = new Coordinate();
                    Coordinate end = new Coordinate();
                    double[] i = vector.GetCoordenatesByCell(1, it);
                    double[] e = vector.GetCoordenatesByCell(vector.Get_area().NumColumns, it);
                    ini.X = i[0] - cellsize;
                    ini.Y = i[1];
                    end.X = e[0] + cellsize;
                    end.Y = e[1];


                    LinearRing line = new LinearRing(new List<Coordinate>() { ini, end });
                    if (test)
                    {
                        IFeature f = li.AddFeature((IBasicGeometry)line);
                        f.DataRow["id"] = it.ToString();
                    }
                    if (fea.Intersects(line))
                    {
                        IFeature rfea = fea.Intersection(line);
                        IList<Coordinate> cor = rfea.Coordinates;
                        if (cor.Count == 1) break;
                        Filldata(rfea, cor);
                        if (test)
                        {
                            IFeature fd = li.AddFeature((IBasicGeometry)rfea);
                            fd.DataRow["id"] = "c";
                        }

                    }

                }


            }
        }

        private void Filldata(IFeature rfea, IList<Coordinate> cor)
        {
            for (int ng = 0; ng < rfea.NumGeometries; ng++)
            {

                Coordinate pi = cor[ng * 2];
                Coordinate pe = cor[ng * 2 + 1];
                int[] cell = vector.GetCell(pi.X, pi.Y);
                int[] cell2 = vector.GetCell(pe.X, pe.Y);

                if (cell[0] == -1 && cell2[0] != -1 && cell[1] != -1)
                {
                    cell[0] = 0;
                }

                if (cell[0] != -1 && cell2[0] == -1 && cell[1] != -1)
                {
                    cell2[0] = vector.Get_area().NumColumns - 1;
                }

                if (cell[0] == -1 && cell2[0] == -1 && cell[1] != -1)
                {
                    cell[0] = 0;
                    cell2[0] = vector.Get_area().NumColumns - 1;
                }

                if (cell[0] == -1 || cell[1] == -1 || cell2[0] == -1)
                { }
                else
                {

                    int row = cell[1];
                    int coli = cell[0];
                    int colf = cell2[0];

                    for (int id = coli; id <= colf; id++)
                    {
                        if (row > 0 && id > 0 && row < vector.Get_area().NumRows && id < vector.Get_area().NumColumns)
                            data[id, row] = 1.0;
                    }
                }
            }
        }
    }

}
